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ABSTRACT 

This  paper  presents  an  approach  to  full-vehicle  simulator 
control  which  accounts  for  nonlinearities  in  a 
vehicle/simulator  system.  The  control  scheme  presented 
is  based  on  the  estimation  of  the  system  inverse 
dynamics.  A  composite  linear/nonlinear  approach  to 
inverse  system  identification  (SYS-ID)  is  presented.  The 
linear  portion  of  the  SYS-ID  uses  time-domain  methods 
to  estimate  the  impulse  response  of  the  inverse  system 
in  a  least  squares  sense.  These  results  are  then 
extended  by  using  the  regularized  approach  to  least 
squares  estimation.  The  nonlinear  part  uses  the  support 
vector  machine  to  approximate  the  nonlinear  deviations 
from  the  linear  model.  Two  approaches  to  using  this 
composite  model  are  presented.  Examples  of  the  linear 
SYS-ID  techniques  are  shown  for  a  2x2  system. 

INTRODUCTION 

With  an  ever  increasing  emphasis  on  vehicle  reliability 
both  in  the  military  and  automotive  industry,  the  science 
of  simulation  and  laboratory  testing  has  correspondingly 
developed.  To  gain  confidence  that  a  particular  vehicle 
will  endure  its  expected  service  environment,  durability 
tests  are  performed  prior  to  production  and  fielding.  The 
time  and  cost  of  such  tests  has  motivated  the 
development  of  laboratory  based  full  vehicle  test  rigs 
(FVTR)  which  have  displaced  much  of  the  time- 
consuming  proving  ground  durability  tests. 

These  FVTRs  were  initially  tire-coupled  and  over  time 
evolved  into  multi-axial  spindle-coupled  configurations 
which  are  able  to  impose  multiple  forces  at  each  spindle. 
Along  with  the  development  of  the  hardware  has  been  a 
corresponding  development  of  the  control  strategies 
used  to  simulate  actual  road  inputs.  The  earliest 
methods  used  road  profile  information,  “effective  road 
profiles”  [10],  or  a  stationary  random  process  to  simulate 
road  roughness.  These  methods  used,  to  some  degree 
or  another,  the  concept  of  the  road  profile  as  a 
conceptual  arbitrator  between  the  test  rig  command  input 
and  the  on-vehicle  response.  Then  Cryer,  Nawrocki  and 
Lund  [5]  removed  the  conceptual  necessity  of  the  road 
profile,  by  modeling  the  relationship  between  the 
command  input  and  the  on-vehicle  response  as  a 


frequency  response  function  (FRF).  They  then  showed 
how  this  FRF  could  be  used  to  estimate  the  proper 
simulator  command,  given  a  desired  on-vehicle 
response.  Their  method  became  the  foundation  of  the 
industry  standard  approach  to  simulator  drive 
determination. 

Although  this  linear  FRF-based  approach  has  been  the 
industry  standard  for  nearly  three  decades,  it  has 
difficulty  compensating  for  nonlinearities  inherently 
present  in  automotive  systems.  To  overcome  system 
nonlinearities,  FRF-based  methods  use  an  iterative 
process  to  converge  on  the  proper  drive  command.  As 
the  techniques  have  been  applied  to  increasingly 
complex  road  simulators,  they  have  been  extended  to 
the  non-square  case  by  Fash,  Goode  and  Brown  [6],  and 
have  incorporated  singular  value  decomposition 
techniques  to  handle  ill-conditioned  FRFs,  but  they  have 
still  remained  dependent  on  the  linear  mathematics  of 
the  FRF  model. 

This  paper  presents  an  approach  to  drive  command 
development  which  incorporates  both  linear  and 
nonlinear  modeling  methods  to  directly  learn  the  inverse 
dynamics  of  the  simulator/vehicle  combination.  First,  it 
describes  the  problem  of  drive  development  in  practical 
and  historical  terms.  Then,  it  discusses  the  mathematics 
and  algorithms  of  the  existing  approaches  and  their 
associated  limitations.  It  formally  states  the  problem  of 
drive  command  development  in  mathematical  terms  and 
gives  the  rationale  behind  the  composite  model.  It  then 
discusses  alternative  approaches  to  both  the  linear  and 
nonlinear  modeling  and  justifies  the  particular  choices 
made.  Finally,  it  presents  some  proposed  alternatives  for 
drive  command  correction.  The  paper  concludes  with 
some  examples. 

BACKGROUND 

The  fundamental  problem  of  full-vehicle  simulation  is  that 
of  replication  of  the  service  environment.  The  service 
environment  is  typically  approximated  by  a  drive  cycle  on 
a  set  of  courses  at  a  proving  ground.  At  the  proving 
ground  vibrational  excitation  comes  from  a  terrain  profile 
represented  by  p(x)  (where  x(t)  is  the  distance  down 
the  course  as  a  function  of  time)  which  creates 


Report  Documentation  Page 

Form  Approved 

OMB  No.  0704-0188 

Public  reporting  burden  for  the  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington 

VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  a  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  OMB  control  number. 

1.  REPORT  DATE 

10  JAN  2005 

2.  REPORT  TYPE 

N/A 

3.  DATES  COVERED 

4.  TITLE  AND  SUBTITLE 

A  Composite  Linear  and  Nonlinear  Approach  to  Full-Vehicle  Simulator 
Control 

5a.  CONTRACT  NUMBER 

5b.  GRANT  NUMBER 

5c.  PROGRAM  ELEMENT  NUMBER 

6.  AUTHOR(S) 

Mark  J.  Brudnak 

5d.  PROJECT  NUMBER 

5e.  TASK  NUMBER 

5f.  WORK  UNIT  NUMBER 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

USA  TACOM  6501  Ell  Mile  Road  Warren,  MI  48397-5000 

8.  PERFORMING  ORGANIZATION 

REPORT  NUMBER 

14118 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

10.  SPONSOR/MONITOR'S  ACRONYM(S) 

TACOM  TARDEC 

11.  SPONSOR/MONITOR'S  REPORT 
NUMBER(S) 

12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release,  distribution  unlimited 

13.  SUPPLEMENTARY  NOTES 

U.S.  Government  Work;  not  copyrighted  in  the  U.S.  Presented  at  SAE  World  Congress  2005,  The  original 
document  contains  color  images. 

14.  ABSTRACT 

15.  SUBJECT  TERMS 

16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION  OF 
ABSTRACT 

SAR 

18.  NUMBER 
OF  PAGES 

12 

19a.  NAME  OF 
RESPONSIBLE  PERSON 

a.  REPORT 

unclassified 

b.  ABSTRACT 

unclassified 

c.  THIS  PAGE 

unclassified 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std  Z39-18 


acceleration  responses  a,  (t)  at  the  spindles  of  the 
vehicle  under  test  (as  shown  in  Figure  1).  On  a  full- 
vehicle  simulator  the  command  inputs  c;(r)  directly  affect 
the  actuator  displacements  df(t) ,  which  in  turn  affect  the 
acceleration  responses  at (t)  at  the  spindles  of  the 

vehicle  under  test  (as  shown  in  Figure  2).  The  standard 
of  validity  for  full  vehicle  simulation  is  that  the  rig 
responses  a,(r)  approximately  match  those  recorded  at 
the  proving  ground  at{t),  which  is  mathematically  stated 
as 

=  Vi  =  l,...,n.  (1) 

Now  the  determination  of  the  simulator  command  c(- (t) 
which  assures  that  this  is  the  case  is  called  drive 
command  development. 

Prior  to  the  work  of  Cryer  et  al.,  drive  file  development 
focused  on  the  measurement  and  replication  of  the 
terrain  profile  p(x) .  The  drive  command  was  then 
derived  from  the  profile  as 

ci(t)  =  fp,i{p(xi(t))) 

for  some  distance  function  xt (t) .  The  function  fp  i( •)  is 

then  used  to  assure  that  (1)  holds.  The  dependence  on 
the  profile  p(x)  is  the  fundamental  limitation  of  this 
approach. 

The  innovation  of  Cryer  et  al.  was  to  directly  replicate 
at{t)  (called  the  desired  response)  by  defining  the  drive 

command  c,(r)  as  being  solely  dependent  on  the  set 

{a,-(f)}"=1  as 

ci (0  =  fc,i  (fli (rt),...,a„ (J,)]  (2) 

where  the  interval  rt  =(t-tup,t-tlow}  denotes  a 
dynamic  relationship  between  the  desired  response 
{ai(t)}"=l  and  Cj(t) .  Let  the  simulator  be  governed  by 
the  dynamics 

«/  0)  =  f a, i  (q  (gt ),  •  •  • ,  c„  ig,))  (3) 

for  gt  —  ^ — oo?  t) .  Equations  (2)  and  (3)  may  then  be 
vectorized  as 


xx{t)  x2(t) 

Figure  1.  Field  environment  to  be  replicated. 


a(0  =  fo(c(?f)) 

(4) 

c(0  =  fc  (a(rt)) 

where  a(t)  =  [dl(t)  ■■■  dn(t)f ,  a(t)  =  [a, (t)  ■■■an(t)]T ,  and 

c(t)  =  [q(f) •••  cn(t)]T  are  nx  1  column  vectors  (i.e.  K'1 ) 

by  convention.  Cryer  et  al.  then  regard  the  system  as  a 
black  box  relationship  between  its  inputs  and  outputs 
without  regard  to  the  unnecessary  concept  of  the  terrain 
profile.  Their  method  will  be  elaborated  in  the 
“EXISTING  METHODS”  section. 

CONTROL  PROBLEM  -  Given  the  definitions  in  (4),  the 
function  fc(-)  must  be  defined  so  that  a (t)  tracks  the 
desired  response  a (t) .  Cast  in  these  terms,  the 
determination  of  the  proper  fc(  )  is  a  control  problem. 

Traditionally  this  control  problem  has  been  solved  by 
attempting  to  guarantee  that 

a(f)  =  a(t)  =  fa  °fc  (a(r)) 
and  solving  for  fc(  )  yields 

W=f aH)-  (5) 

If  f„(-)  is  perfectly  known  and  is  invertible  everywhere 
over  its  range  then  (5)  assures  that  (1)  will  hold.  Now 
given  that  fa(  )  consists  of  the  dynamics  of  vehicle 
specimen,  the  simulator,  and  the  PID  controller,  it  is  not 
readily  derivable  from  first  principles.  Additionally,  fa  (■) 

is  typically  not  known  down  a  finite  number  of  unknown 
physical  parameters  due  to  the  presence  of  the 
unmodeled  vehicle  dynamics.  Worse  still  fa(-)  is 

nonlinear.  Given  these  circumstances,  existing  practice 
is  to  perform  a  system  identification  on  the  plant  (vehicle 
and  simulator)  to  obtain  the  estimated  system  dynamics 
ffl(-) .  The  following  section  briefly  describes  the  method 
of  Cryer  et  al.  in  the  determination  of  ft.  (•) . 

EXISTING  METHODS  -  The  methods  developed  by 
Cryer  et  al.  are  based  on  discrete-time,  linear  math  in  the 
frequency-domain.  Let  the  input  and  output  signals 


Figure  2.  Simulator  with  HMMWV  attached. 


(7) 


(q(t)  and  ^(0  respectively)  be  discretized  at  a  sample 
time  Tj  such  that  t  =  k ts  .  Also  let  the  inputs  and 
outputs  be  transformed  into  the  frequency-domain  as 
aifj)  =  S{a(k)}  and  c(fj)  =  ${c{kj)  where  ${•} 

denotes  the  discrete  Fourier  transform  (DFT).  Cryer  et 
al.  then  define  fa  (•)  as 


c,+1  (k)  =  c,  (k)  +  wi+1  o  (Hinv  (•)  *  e(-  (*)) 


and  progress  is  measured  by  monitoring  a  statistic  of  the 
error  which  is  typically 


RMSt(e,-0) 
RMS*  (a(-)) 


(8) 


a(*)  =  fo(c(fc))  =  H0*c(*)  (6) 

where  H(&)  is  the  system  impulse  response  function 

and  *  denotes  the  convolution  operation.  Equation  (6) 
has  the  corresponding  frequency-domain  representation 

a(fj)  =  H(fj)c(fj) 

where  H (/■)  is  the  system  FRF.  Since  H (/•)  is  not 

known,  it  must  be  experimentally  determined,  so  let  n(k) 
be  a  suitable  colored-noise  signal  used  to  excite  the 
simulator,  obtaining  the  response  r (k) .  If  n (k)  and  r(k) 

have  DFTs  n(/;)  and  r(/;)  respectively,  then  the 
frequency  response  function  is  estimated  as 

Wj)  =  Srn(fj)Snl(fj) 

where  Snn(/;)  =  E  n(fj)nH(fj)  is  the  auto-spectral 

density  (ASD)  of  n (k)  and  Srn(fj)  =  E  r(fj)nH(fj)  is 
the  cross-spectral  density  (CSD)  of  r (k)  and  n(k) .  Now 
according  to  (5),  H (/•)  must  be  inverted  to  obtain  the 
estimated  drive  c (k)  from  a (k) ,  thus  yielding 


where  RMS^.  (•) :  RnXt  R"  and  division  is  computed 
element-wise.  This  iterative  process  is  terminated  when 
each  element  of  m(  is  sufficiently  small. 

In  practice  this  iterative  process  requires  10  to  20  such 
iterations  to  obtain  satisfactory  convergence.  It  is  widely 
known  that  the  success  of  the  iteration  process  is 
dependent  on 

1.  The  quality  of  the  forward  model  H(&). 

2.  The  choice  of  excitation  n (k) . 

3.  The  choice  of  the  weights  w(- . 

4.  The  stationarity  of  the  desired  response  a (k) . 

Each  of  these  issues  can  cause  difficulty  in  obtaining 
sufficiently  small  relative  errors  m,  and  make  the 
process  highly  dependent  on  the  skills  of  the  practitioner. 
It  is  believed  that  most  of  these  difficulties  originate  from 
the  essential  nonlinear  nature  of  the  plant.  Thus  to 
obtain  effective  improvement  in  drive  file  development, 
we  claim  that  the  nonlinearities  of  the  plant  must  be 
incorporated  into  the  inverse  model  fc(  ).  It  is  this  goal 
that  motivates  the  approach  presented  in  the  sequel. 

PROBLEM  DEFINITION  AND  APPROACH 


c(0  =  fc(a(0)  =  Hinv0*aCfc) 

where  Hinv(0  =  ^_i{h_1(/7)}  .  In  the  frequency- 
domain  the  drive  command  estimate  is 

c(fj)  =  K-\fj)a(fj) 

which  is  independent  for  each  frequency  line  fj  . 

Now  because  H(£)  is  a  linear  estimate  of  a  nonlinear 
system,  c (k)  will  not  be  optimal.  To  obtain  the  best  drive 
estimate  the  iteration  process  is  used.  Let  the  initial 
drive  estimate  be  defined  as  c0(£)  =  w0  °(Hinv  0*a(*)) 

where  w0e[0,l]"  is  a  weight  vector  and  °  denotes  an 
element-wise  product.  Then  c0(k)  is  used  to  drive  the 
simulator  and  the  corresponding  response  a0(£)  is 
recorded.  Let  e0(k)  =  a(k)-a0(k)  denote  error  between 
the  desired  response  and  actual  response,  then  the  drive 
estimate  correction  is  given  by 

Ci(*)  =  co(*)  +  w1°(Hmv0*eo(*)) .  This  process  is 
repeated  according  to  the  rule 


Since  the  forward  dynamics  of  the  plant  fa(  )  are 
nonlinear,  the  inverse  fc(  )  will  also  be  nonlinear. 

However  since  the  linear  model  Hinv(0  is  used  to  yield 
a  first  order  approximation  of  the  inverse,  it  is  clear  that  it 
contains  significant  information  about  the  plant 
dynamics.  Therefore,  the  linear  model  will  be  retained 
here  and  a  nonlinear  modeling  approach  will  be  used  to 
approximate  the  deviations  from  the  linear  model.  This 
approach  is  illustrated  in  Figure  3  where  the  forward 
dynamics  are  modeled  as 

fa(c(0)  =  H0*c(0  +  (a(c(^))  (9) 

and  the  inverse  dynamics  are  modeled  by 

f,  (a(*))  =  Hinv  (•)  *  a (k)  +  V  («(** ))  (10) 


Forward  Model  Inverse  Model 


Figure  3.  Linear/Nonlinear  approach  to  modeling. 


(13) 


the  set  of  ordered  pairs  {(;',  y)}"’^  <->A.  This  “flattened” 
version  then  takes  the  form 


jjinv  a 


where  h •/nv7’ =  r^i/nv (fcmin )  •••  /z;ylv(£max)J  is  the  inverse 
impulse  response  between  dj(k)  and  clin.(&).  Let 
Rnr  D<pa(£)  =  {a (k-Tj)}k™?  which  is  indexed  by  A  and 

—  V  "-min 

r  ~\T 

given  by  <pa(£)  =  \aT(k-kmin)  •••  a T(k-kmax)  . 

Equation  (11)  may  then  be  expressed  as  a  matrix 
product 


clin(£)  =  Hin>a(£)  (12) 


and  will  be  used  as  a  basis  for  estimating  Hinv(£)  in  the 
following  section. 


LEAST  SQUARES  METHOD  -  Now  as  stated  earlier, 
the  system  is  excited  with  colored  noise  n(k)  and  the 
response  r (k)  for  k=\,...,t  is  recorded.  Then  for  each 
r (k)  we  form  the  regression  vector  <pr  (A)  using  the 

index  A  where  clearly  r(i)  =  0  for  k<\  and  k>l.  By 
hypothesis,  the  inverse  plant  model  obeys  (12)  for  any 
set  of  inputs  and  outputs  which  yields  the  transposed 
relation 

n  ( k )  =  (pr  (A:)H 


Let 

N  =  [n(l)  •••  nCQfeM601 

and 

®r=[?r(l) 

~\T 

■■■  <pr  (O  eRexnr;  then  the 

following 

aggregated  model  is  obtained 

N  =  <DrHinvr  . 

The  Moore-Penrose  pseudo-inverse  may  be  used  to 


solve  for  Hinv  as 

Hinvr=(®^®r)  'oJn 

which  may  be  shown  to  minimize  the  summed  loss 

J (Hmv )  =  ^X  “(*) - Hm>r (*)  2  (14) 

2k= l  2 

and  hopefully  the  expected  loss 

2" 

E  n(£)-Hinv<pr(T)  .  (15) 

REGULARIZED  LEAST  SQUARES  METHOD  -  Now 
Hinv  as  learned  using  the  least  squares  method  is  a 

non-parametric  estimate  with  n2r  free  parameters. 
Generally  the  half-length  of  the  impulse  response 
max ( | Amin  | , |^max  |)  is  chosen  to  correspond  roughly  to  the 

settling  time  of  the  system.  So  for  a  small  sample  time 
Tj  and  long  settling  time  we  typically  have  kmax  ~  100 

which  means  that  Hmv  will  have  hundreds  if  not 
thousands  of  free  parameters.  This  is  typically  more 
than  the  number  of  degrees  of  freedom  nd  in  the  actual 

system.  Therefore  an  ideal  model  would  only  contain  nd 
degrees  of  freedom.  Unfortunately,  nd  is  not  known  but 

it  may  be  asserted  with  confidence  that  n2r  »  nd .  Due 
to  this  over  parameterization  and  that  (13)  regards  the 
elements  of  Hinv  as  being  independent,  the 
phenomenon  of  over  fitting  is  expected  which  means  that 
the  minimization  of  (14)  does  not  generally  imply  good 
performance,  i.e.  small  (15).  Such  a  problem  with  more 
parameters  than  degrees  of  freedom  is  called  ill-posed. 

The  solution  of  ill-posed  problems  was  first  studied  by 
Tikhonov  et  al.  [8]  who  found  that  the  phenomenon  of 
over-fitting  can  be  mitigated  by  introducing  a  so-called 

regularization  functional  -p(Hinv) :  Rnxnr  M+  into  the 

optimization  problem  (14)  which  will  penalize  those  Hinv 
which  are  least  preferred.  Now  it  is  required  that  V(-)  be 


Figure  5.  Forward  Impulse  Response.  Figure  6.  Inverse  Impulse  Response. 


positive  semi-definite  and  is  a  super  set  of  the  null  space 
of  J(-)  in  (14).  With  these  two  conditions  met,  it  may  be 
chosen  to  convey  a  preference  for  a  particular  structure 
of  Hinv .  Thus  the  minimization  problem  now  has  two 
objectives  which  must  be  mutually  weighted  with  the 
parameter  y  as 

minimize:  yV  ( Hinv )  +  J  ( Hinv ) . 

jjinv  \  /  \  / 

Commonly  the  regularization  functional  is  chosen  to  be 
of  the  following  quadratic  form 

Vq  (Hinv )  =  |  vec  (Hinv  j7  Q  vec  (Hinv )  (16) 


where  vec(A) :  Rmxn  Rmn  vectorizes  its  argument  for 

some  AeMmx"  by  stacking  its  columns  and 

2  2 

Q  =  Qr  e  R"  rxn  r  is  a  positive  semi-definite  symmetric 
matrix.  The  regularized  solution  then  becomes 


H 


invr 


-1 


(l„®®r®r  +  TQ)  1  ®  ®r  )  vec(N)  j 


(17) 


where  ®  denotes  the  Kronecker  product  of  two  matrices 
and  vec-1(-)  reverses  the  effects  of  vec(-).  Now  we  note 
if  Q  =  I,  then  (17)  is  known  as  ridge  regression  and 
expresses  a  preference  for  a  minimum  norm  solution.  In 
the  case  of  time-domain  modeling  there  are  other  more 
useful  choices  for  Q  as  will  be  developed  next. 


Time-Domain  Regularization  -  As  mentioned,  the  matrix 
Q  may  be  chosen  to  express  a  preference/penalty  for 
certain  types  of  impulse  responses.  This  section 
describes  how  to  choose  Q  to  express  preferred  time- 

domain  behavior  of  Hinv.  Without  loss  of  generality 
consider  the  SISO  (i.e.  n  =  1)  system  with  inverse 
impulse  response  him{k)  and  where  Hinvr  =  hinv  e 


is  a  column  vector 


hmv(kmi  n) 


hmV(kmax) 


~\T 


which 


extends  forward  and  backward  in  time.  Due  to  this 
acausality,  it  is  expected  that  lim  hmv(k)  =  0.  If  kmax 

k— >±oo 


roughly  corresponds  to  the  settling  time  of  the  system, 
then  for  exponential  decay  it  is  reasonable  to  expect  that 


hmv(k) 


-k 


<  e{k)  -  or  exp  -TJ  for  some  a  and  kc  (see 


Figure  7).  Then  if  the  expected  settling  behavior  is  such 


that 


i  inv  /  / 

h  (kn 


h'nv  (0) 


where  0<//<sl  then  a 


regularization  functional  Vk  (hinv)  =  ihinv7hmv  would 
unfairly  penalize  those  values  for  which  k~  0  and  will 
tend  to  equalize  all  values  of  hmv(k).  If  instead  the 


values  of  /zinv (A)  were  scaled  with  the  envelope  e(k) 
then  the  settling  preference  would  be  properly 
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Figure  7.  Decaying  Impulse  response. 


expressed.  The  regularization  functional  is  then  defined 
as 


n 


,  .  \  1  ^max 

(/zinv(-)) = —  £ 

'  iinw  tv 

\  /  9  ^ 

^  k—k  ■ 
K~Kmm 

{  e(k)  j 

-  hinvTQ;.hil 


(18) 


where  [QJ..=-^ — — - =  &7expf2^  A’min+1l\  Then  by 

L  e2(i-kmia+ 1)  11  kc  ) 

choosing  kc  (which  can  be  thought  of  as  a  time 
constant),  the  preferred  settling  behavior  of  h(k)  may  be 
expressed  via  Vk  (•) . 


Frequency-Domain  Regularization  -  Similarly, 
preferences  for  the  frequency-domain  behavior  of  h(k) 

(we  drop  the  super-script  (-)inv  for  ease  of  notation)  may 
be  expressed  using  a  regularization  functional.  First 
observe  that  h(k)  is  a  discrete-time  aperiodic  signal 
which  uses  the  following  Fourier  transform  pair 


H(e)=S{h(kj}=  Y,  h(k)e~'iek 

k=—oo 
1  71 

h(k)  =  3  '{H(0)}  =  —  J  H(d)ei6kdd 

-71 

where  j  =  V4  and  de[-K,rr)  corresponds  to  angular 
frequency  such  that  an  unsampled  frequency  may  be 
obtained  as  /  =  ^TS  where  fv  =  J- .  To  penalize  h(k) 

in  the  frequency-domain,  let  W(6) :  [-rr.rr)  M+  be  a 

positive  symmetric  weighting  function  which  ascribes  a 
penalty  to  particular  values  of  H{9)  for  every  6.  A 

frequency-domain  regularization  functional  may  then  be 
defined  as 


Vd{H (•))  =  ~  j  \W(9)H(9)\2  dd  (1 9) 

-71 

which  is  expressed  in  terms  of  H(0)  =  ${h(k)}  which  is 

not  typically  available.  Using  Parseval’s  Theorem  (19) 
may  be  converted  to  the  time-domain  as 


Mk-)) 


2 


i=k„ 


^max  00 

2  h(i)h(j)  w(k)w(k  +  i-j ) 


J=kn 


k=—oo 


where  h’(^)  =  g?'1{W(6')}  has  infinite  extent  (i.e. 
k  g  ( — °°, 00  ) )  but  does  not  have  to  be  calculated  explicitly. 

Let  [Q/]..  ='£J™=_00w(k)w(k  +  i -  j) ,  then  it  may  be 

shown  that  ^  Z-00  Mk)w(k  +  i  -  j )  =  3-'  {w2  (0)}  (. i-j ) 

then 


Ve  (AO)  =  {hrQ0h 

where  Q0  is  symmetric  and  positive  definite  and  may  be 

determined  if  $'lw2(0)]  can  be  calculated.  A  typical 

example  of  this  technique  is  when  h(k)  should  be  limited 
to  the  band  of  frequencies  9km  <  0  <  .  In  this  case 

W(6 >)  =  • 


up  ■ 

0  0mv<e<e„ 

1 


7low  —  v  —  ^up  ’ 

otherwise 


then  it  is  easily  shown  that 

sin  ( 9up  (i-j))- sin  ( 6\ow  (i  -  j) ) 


[Q<?L  =  sa  ~j)~- 


for  i*j  and  [QJ..  =1- 


Aip  ^low 


Mi- j ) 
for  i  =  j  ■ 


Combined  Regularization  -  These  two  regularization 
methods  are  often  used  together.  They  may  then  be 
combined  as 


Q  -  aQk  +^Q<? 

where  if  a,b>  0  then  Q  is  then  positive  semi-definite. 


NONLINEAR  MODELING 


CONSIDERATIONS  -  Given  the  desire  to  find  a 
nonparametric  model  \j/C)  a  suitable  approach  must  be 
chosen.  Several  techniques  are  available  to  estimate 
i|/(-) .  The  final  selection  will  be  based  on  the  following 
considerations. 


Capacity  Control  -  When  one  seeks  to  learn  a  general 
mapping  i|/(-) ,  all  nonparametric  methods  use  a 
particular  fixed  structure  which  is  parameterized  by 
kl5  as  \|/()  =  »)/(•;  k)  such  that  the  following 
optimization  problem  is  solved 

minimize:  J(^)  =  -XL(nnonUn(^X'l'(?r(^);^))  (21) 

1  L  k= 1 

given  a  loss  function  L(-,-):R"xR''hR+  which 

generates  a  penalty  for  \|/(<pr  (£);}.)  *  nnonlin  (£) .  It  is 

widely  known  that  the  ability  to  minimize  (21)  over  a 
training  set  (20)  is  dependent  on  the  number  of  free 
parameters  k,  the  structure  of  v| /(•;■).  the  inherent 

complexity  of  the  process  which  generated  the  training 
data  (20),  and  the  size  of  the  training  set  l .  Given  these 
concerns  it  is  important  to  choose  v|/(-;X)  to  avoid  the 
phenomenon  of  over  fitting  defined  as 


-  “nonlinW'V(?rW;^))  »  ^  (K) 


where  k0  is  the  optimum  obtained  in  (21).  Vladimir 
Vapnik  [9]  formalized  the  notion  of  “freedom  to  fit”  by 
defining  the  notion  of  capacity  of  X.)  to  fit  a  set  of 

data.  This  capacity  is  theoretically  quantifiable  using  the 
VC-dimension,  A  .  The  phenomenon  of  over-fitting  is 
then  to  be  expected  for  n>£.  It  is  desirable  then  to 
have  an  estimator  with  a  low  VC-dimension.  Generally 
one  should  expect  reasonable  performance  for  an 
estimator  for  which  t  >  20n  . 


Once  the  linear  portion  of  the  inverse  model  is  learned, 
then  the  residues  become  the  nonlinear  part  of  the 
inverse,  therefore 

nnonlin(k)  =  n(k)-Hin\r(k) 

and  the  nonlinear  modeling  task  consists  of  learning  a 
relationship  vy :  K,,r  k>  R"  such  that 

\j/  (<pr  (A))  =  nnonlin  f A)  +  e(A) .  Because  a  nonparametric 

technique  was  used  for  the  linear  portion  of  the  model,  it 
is  appropriate  to  use  a  nonparametric  technique  for  the 
nonlinear  portion  as  well.  The  data  set 

{fer^XnnonlinW)}^  (20) 

is  available  with  which  to  train  the  estimator. 


The  “Curse  of  Dimensionality”  -  Due  to  the  nature  of 
nonlinear  relationships,  many  more  training  samples  are 
required  to  learn  a  nonlinear  relationship  than  a  linear 

one.  In  general,  to  learn  a  mapping  it  is 

necessary  to  observe  its  behavior  in  all  “regions”  of  the 

bounded  domain  DcR^.  The  domain  V  must  then  be 
decomposed  into  a  set  of  regions  Q  in  which  its 
behavior  must  be  observed.  To  properly  estimate  the 
mapping,  each  Q,  must  be  sufficiently  small  such  that  a 

sample  xe  Q  contained  therein  represents  all  of  L2,- 
and  the  set  U,T2(-  must  cover  V  .  With  a  vector-valued 
domain,  the  number  of  regions  necessary  to  decompose 
the  domain  is  exponential  in  the  dimension  d  of  the 
domain.  So  as  the  dimension  of  the  input  space 
increases,  it  becomes  increasingly  difficult  to  obtain 
sufficient  samples  to  fill  V .  This  is  called  the  “curse  of 
dimensionality”.  The  reader  is  referred  to  Figure  8.  This 
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Figure  8.  The  curse  of  dimensionality. 


Rd 


reality  is  a  fundamental  limitation  of  non-linear 
approaches.  It  motivates  the  user  to  limit  the  dimension 
of  their  input  space  as  much  as  possible  and  limit  the 
extent  of  the  domain  to  the  smallest  possible  set. 


where  |-|  is  called  the  e  -insensitive  loss  function  which 

is  defined  as  \e\  =  max(0,|e|-f)  and  is  illustrated  in 
Figure  9.  The  optimum  estimator  (22)  is  then  given  by 


-e  +s  e 

Figure  9.  e  -insensitive  loss  function. 


Available  Approaches  -  Given  the  problem  of  learning  a 
general  mapping  many  approaches  are  available. 
Methods  which  have  been  studied  in  the  literature 
include  artificial  neural  networks  [1],  Takagi-Sugeno 
fuzzy  logic,  wavelets,  splines,  etc.  Recently  a  new  class 
of  methods  called  kernel  methods  has  been  shown  to 
have  some  advantages  over  these  other  methods.  By 
far  the  most  prominent  kernel  methods  to  date  are  the 
support  vector  machines  (SVM).  Their  desirable 
features  include: 

1 .  They  generalize  well. 

2.  They  are  based  on  linear  mathematics. 

3.  They  result  in  a  global  optimum. 

4.  They  are  solved  using  quadratic  programming 
techniques. 

The  support  vector  machine  will  be  used  to  learn  the 
mapping  y :  R"r  R"  on  the  set  given  in  (20). 

SUPPORT  VECTOR  MACHINE  -  The  support  vector 
machine  is  a  recently  developed  method  for  solving 
nonlinear  classification  and  regression  problems  (see 
Cristianini  et  al.  [4]).  It  was  invented  by  Vapnik  and 
extended  by  many  others.  The  following  two  sections 
present  a  synopsis  of  the  SVM  regression  problem  for 
the  scalar-valued  and  vector-valued  cases. 

Scalar-Valued  Support  Vector  Regression  -  Given  a 
general  mapping  /:  R'"  R  and  the  training  set 

y-t )}._!  the  SVM  seeks  to  find  an  estimator 

y{  =  f(Xj)  which  will  perform  well  on  the  training  set  and 

is  expected  to  perform  well  in  use.  The  SVM  estimator 
takes  the  following  form 

y(x)  =  (w,<p(x)}  +  fc  (22) 

which  is  linear  in  the  mapped  input  vector  <p(x,)  and  has 
free  parameters  w  (the  weights)  and  b  (the  bias).  The 
mapping  <p :  R”  k>  Rv  takes  any  vector  from  the  input 

space  R”  to  a  higher-dimensional  feature  space  Rv , 
that  is  v»n.  To  determine  the  parameters  w  and  b 
which  are  optimal,  a  loss  function  is  necessary.  In  the 
standard  SVM  formulation, 


the  weight  and  bias  which  minimizes 
J(w,b)  =  |y,-  —  y,-]  -  It  is  expected  that  v>£ 

implying  that  the  optimum  w  and  b  are  not  unique. 
Vapnik  has  shown  that  small  ||w||  usually  generalize 

better  than  large  ||w||  (that  is  “flatter”  models  generally 
predict  better).  As  such,  the  SVM  method  employs  the 
regularization  functional  V(yr)  =  4||wf  yielding  the  so 
called  primal  optimization  problem 

0  minimize:  i||w| |2 +C^=1|yf 

where  %  is  defined  in  (22).  Because  of  the  potentially 
high  dimension  v  of  the  feature  space,  the  problem  is 
recast  into  dual  form  (see  Bazaraa  et  al.  [2])  as 

i  i  i 

\d\  maximize :  -  \  Z  Mjk^i  ’ x; )  +  Z  >'<A  “  ^Z I A I 

{#},=!  U=i  <=i  i=i  ^3) 

such  that  XA=0  ancl  \Pi  \  -  C 
(=1 

where  k(xhxj)  =  ^<p(x,-),q>(x;-)^  is  the  so-called  kernel 

function  which  defines  an  inner  product  in  feature  space. 
The  weight  vector  in  the  dual  variables  becomes 

w  =  Z;=1A<P(xi)  anc^  f°rm  estimator  then 

becomes, 

t 

KX)  =  Z PM^i,s)  +  b  ■  (24) 

(=1 

Notice  that  both  the  optimization  problem  (23)  and  the 
estimator  (24)  are  expressed  in  terms  of  the  kernel 
function  k(-,-)  only,  not  the  mapping  <p() .  Mercer  [7] 

showed  that  a  function  which  satisfies 

J  g(x)k(x,y)g(y)dxdy  >0,  Vg(-)e  L2(T) 

XxX 

for  X  =  dom(^())  implicitly  defines  the  mapping  <p(  ) , 

and  guarantees  that  it  is  an  inner-product  space.  Such  a 
kernel  is  called  a  Mercer  kernel. 


The  solution  to  (23)  is  typically  sparse  so  let 
Is  ={/:/?,  =£  0}  be  the  index  set  of  non-zero  dual 
variables.  Then  the  so-called  support  vectors  are 
defined  by  {x(}.gJ  .  This  sparseness  property  of  the 

SVM  allows  them  to  adjust  their  capacity  to  the 
complexity  of  the  training  data.  This  allows  them  to 
generalize  well . 

Vector-Valued  Support  Vector  Regression  -  A 
generalization  of  scalar-valued  SVR  is  the  vector-valued 
SVR  which  is  designed  to  learn  a  mapping  of  the  form 

given  a  training  set  {(x,-,y(- )}f—1  ■  In  VV- 

SVR  (see  Brudnak  [3]  for  full  development)  the  estimator 
takes  the  general  form 

y(x)  =  W<p(x)  +  b 

where  WeRrfxi/  is  a  weight  matrix  and  bel^  is  the 
bias  vector.  The  vector-valued  loss  function  is  defined 
as 

T(y,,y,)=  |y,--y,-|l 

F  £ 

where  1  <  P  <  oo  represents  the  y>-norm  of  the  error. 
The  primal  optimization  problem  is 

0  minimize:  ±||W||*  +^-^1  -y«-||p|£ 

which  leads  to  the  dual  problem 

t  t  i 

[ d\  maximize:  £  r,T;-A-(x,-,x;-)  +  £y,T,- 

{r/},=l  »,/=!  i= 1  i= 1 

l 

such  that  Yr,  =0  and  ||r,  ||  <C 

V—i  i  II  i\\q 

i—l 

where  ||-||  and  IHI  are  conjugate  norms  (i.e.  — +  -  =  1). 

II  Up  w\\q  J  &  pq 

Given  that  W  =  2J  .=iri-<p(xI-) ,  the  final  form  of  the 
estimator  becomes 

t 

y(x)  =  0r,T(x,-,x)  +  b  . 
i= l 

The  index  set  for  the  support  vectors  becomes 
Js  =  j/' :  ||rf||  ±  oj  and  the  support  vectors  are  defined  as 

{xJfej  ■  This  vector  valued  version  is  used  to  learn 
nonlinear  mappings  for  MIMO  systems. 

COMPOSITE  MODELING 

Now  two  models  exist  of  the  inverse  of  the  plant.  One  is 
linear;  the  other  is  nonlinear.  A  composite  model  is 
created  as  shown  in  Figure  3  and  Figure  4  by  summing 
the  output  of  the  two  models.  To  estimate  these  models 
the  plant  is  excited  with  the  random  time  series  n(k)  and 


the  associated  response  r (k)  is  recorded.  From  these, 

the  data  set  {(nffcXrCfc))}^  is  used  to  learn  the  inverse 

impulse  response  function  Hinv(-)  as  discussed  above. 
The  best  linear  estimate  of  the  inverse  plant  is  therefore 
defined  as  Hinv(-).  The  predicted  input  based  on  the 

linear  model  then  becomes  n(t)  =  Hmv()!»r(t).  The 
deviations  from  the  linear  model  may  then  be  estimated 
as  d(k)  =  n(k)-n(k) .  By  hypothesis  d(l)  contains 
information  regarding  the  nonlinear  part  of  the  model, 

therefore  the  training  set  D0  ={(<pr(£),d(&))}^  ^  may  be 

learned  as  a  mapping.  The  vector-valued  support  vector 
machine  approach  will  be  used  to  accomplish  this.  The 

estimated  mapping  will  be  denoted  as  \j>0  :  R"r  R" 

and  takes  the  form  »j/0  |{r(fc)}*;™1  J  =  \j>0  (<pr(£)j .  Once 

this  mapping  is  learned  it  may  be  used  to  estimate  the 
proper  drive  which  will  yield  the  proper  response  as 
follows 

c(&)  =  Hinv9a  (A)  +  Vj>0  (<pa(£)j 

where  q \(k)  is  the  vectorized  (i.e.  indexed  by  A) 

desired  response  a (k).  So  if  Hmvand  \j>0(  )  properly 
represent  the  inverse  dynamics  then  c (k)  is  will  be  a 

good  estimate  of  the  drive  and  it  is  expected  that  the 
actual  rig  response  will  be  equal  to  the  desired  response, 
that  is 

fa(c(*))»a(*). 

If  the  desired  response  a (k)  is  statistically  dissimilar 
from  the  SYS-ID  excitation  response  r (k)  then  the  curse 
of  dimensionality  indicates  that  Vj>0(-)  will  not  necessarily 
be  a  good  model  for  the  inputs  a (k) .  In  such  cases  the 
nonlinear  mapping  may  have  to  be  re-learned  for  plant 
outputs  which  are  closer  to  the  desired  response.  The 
next  section  discusses  a  couple  of  alternative 
approaches. 

THE  CONTROL  SCHEME 

In  this  section  strategies  for  determining  the  best 
estimated  drive  c (k)  are  discussed.  Typically  the  first 

estimate  of  the  drive  (called  c 0(k))  is  not  the  best 
possible  drive  by  the  criteria  stated  in  (8).  To  achieve  the 
optimum  drive,  typically  a  sequence  of  drives  is 
estimated  c0(-), q cf(-)  until  an  acceptable  level  of 

accuracy  is  attained.  Two  alternative  strategies  are 
presented  here. 

PREDICT-CORRECT  -  This  method  requires  a 
linearization  of  the  estimated  inverse  plant.  Let 


then 


fi inv  (9  a  (*))  ±  Hin>a  (k)  +  vo  (<pa  (*)) 

*"[*)*£r  .(#.) 

where  Jmv  (<pa)  is  the  Jacobian  of  f inv  (<pa ) .  The  rig 
response  which  is  generated  by  the  drive  £,-(■)  is  a,-(0 
and  the  associated  error  is  given  by  ef  (■) .  The  Jacobian 
may  then  be  used  to  estimate  a  new  drive  given  by 

«/+ 1  (')  =  c«  (■)  +  wi+1  o  ( Jinv  (<pa )  e,-  (■)) 


2.  If  m,  = - V,  .  '  <  t  for  some  threshold  r 

'  RMS*(a(*)) 

then  stop. 

3.  Let  di(k)  =  ci(k)-Hiny^(k). 

4.  Form  2)i+1  ±  S,  u{(?g.  (fc),d;(*))}^=1  . 

5.  Train  \\ii+l( ■)  such  that  \\ii+ j(-)  i=  Di+1 . 

6.  Estimate  the  next  drive 

ci+i(^)  =  Hinvfa«)  +  yi+1(<paW). 

7.  Play  c i+1(k)  into  the  rig  and  record  the  response 

a«+i  (  k  )  ■ 

8.  Goto  step  2. 


where  w,+1e[0,l]"  is  a  weight  vector.  Because  of  the 

nature  of  the  method,  it  is  called  predict-correct.  First, 
note  that  this  method  reduces  to  the  linear  method  in 
Equation  (7)  when  \j>0(cpa)  =  0.  Secondly,  this  method  is 

expected  to  work  well  if  vj>0(cpa)  is  a  very  good 

approximation  of  the  nonlinear  portion  of  the  inverse 
model’s  derivative.  It  is  not  expected  that  this  will  be  the 
case.  For  this  reason  the  approach  presented  in  the 
following  section  is  preferred. 

OBSERVE-ADAPT  -  Another  method  which  does  not 
require  the  calculation  or  estimation  of  a  derivative  is  the 
method  called  observe-adapt.  In  this  method,  the  linear 

part  Hinv  is  held  fixed  and  the  nonlinear  part  of  the 
model  is  adapted  to  new  plant  observations.  Given  the 
initial  nonlinear  inverse  model  Vj>0(<pa),  which  was 

trained  on  £>0 ,  it  is  used  to  calculate  an  initial  drive 
estimate  c 0(k)  which  produces  an  associated  rig 
response  a0(-) .  It  is  expected  that  this  response  will  not 
be  optimal,  however,  by  hypothesis,  it  should  be  closer  to 
desired  response  a (k)  than  to  the  original  random 

excitation  response  of  r (k) .  In  this  case  then  the 
associated  deviation  from  the  linear  model  may  be 
estimated  as  d0(k)  =  c0(£)-Hinva0(k)  and  if  the  desired 
data  is  indexed  by  (1,...,L)  an  augmented  training  set 
may  then  be  developed  as 

®i  =  (*),d0(*))}^  .  The  set  S),  may  then  be 

used  to  estimate  a  new  mapping  such  that  viq  (■)  1= 
where  the  symbol  h  denotes  that  ^iO)  models  or 
generalizes  the  set  £>j .  This  new  nonlinear  model  may 
then  be  used  to  produce  a  better  drive  as 
£](£)  =  Hinv<pa(£)  +  \j/,  (<pa(T)) .  This  process  may  then 
be  generalized  to  the  following  algorithm. 

1 .  Given:  The  estimated  drive  c {(k) ,  the  corresponding 
response  cp^.  (k),  the  linear  model  Hinv  and  the 
training  set  £>,■ . 


In  this  way,  if  the  sequence  of  responses  continues  to 
get  closer  to  the  desired  response  (i.e.  m;  >mi+1)  then 

the  sequence  of  models  (v0(’)»  ViO).  •••»  V/C).  ■••)  should 

converge  to  the  best  representation  of  the  nonlinear 
behavior  at  the  operating  point. 

EXAMPLE 

To  demonstrate  some  of  the  methods  presented, 
consider  a  2x2  system  (i.e.  n  =  2)  which  comprises  the 
rear  two  tires  of  a  single-axle  military  trailer  (see  Figure 
10)  where  the  two  inputs  cx{k)  and  c2(k)  are 
displacement  commands  to  the  left  and  right  actuators 
respectively  and  the  responses  ax(k)  and  a2(k)  are  the 

associated  acceleration  responses  at  the  left  and  right 
wheel  spindles.  The  control  problem  is  then  to  develop 
the  proper  drive  command  c  (k)  such  that  the 
acceleration  responses  match  some  a-priori  specified 
response  a (k) .  In  this  example  the  linear  SYS-ID 
portion  of  the  above  developed  process  will  be 
demonstrated. 

To  begin  the  system  identification  process  a  random 
noise  command  n (k)  is  generated  with  which  to  excite 
the  system.  Such  excitation  commands  are  typically 
specified  in  the  frequency-domain  as  shown  in  Figure  11. 
They  are  then  converted  to  the  time-domain  such  that 
they  are  uncorrelated  in  the  time-domain.  Because  the 
system  is  not  linear,  it  is  desirable  to  excite  the  system  at 


Figure  1 0.  Military  trailer  which  comprises  2x2  system. 


all  energy  levels  from  mild  to  severe,  the  time  history  is 
linearly  increased  in  amplitude  throughout  its  duration.  In 
this  particular  case  it  was  designed  to  begin  at  roughly 
±0.5  cm  and  progressively  increase  to  roughly  ±2.5  cm. 
The  generated  time  histories  are  shown  in  Figure  12. 
This  time  history  was  played  into  the  simulator  and  the 
associated  response  r (k)  was  recorded.  Both  of  the 
input  and  output  were  discretized  at  a  sample  rate  of 
fs  =  204.8  Hz  The  techniques  of  linear  time-domain 
modeling  were  applied  to  the  recorded  inputs  and 
outputs  to  learn  the  inverse  time-domain  model  Hinv(£) . 
For  reference  purposes,  the  inverse  impulse  response 
function  was  calculated  for  kmiB  =-400  and  kmax  =400 
which  is  shown  in  Figure  13  where  the  average  mean 
squared  error  was  0.1209  cm. 

The  time-domain  regularization  technique  is 
demonstrated  by  training  on  the  same  information 

however  expressing  a  preference  that  Hinv(£)  decay  to 
a  value  of  10%  at  £s=400.  The  estimated  impulse 
response  function  is  shown  in  Figure  14.  In  this  case  the 
average  mean  squared  error  was  0.1218  cm  which  is 
very  close  to  the  unregularized  optimum.  It  is  therefore 
observed  that  the  performance  has  degraded  by  about 
0.7%  but  the  model  is  much  different.  This  is  evidence 
that  the  original  solution  in  Figure  13  was  unstable  or  the 
problem  was  ill-posed. 


Figure  12.  Plot  of  random  excitation  n (k) . 


The  frequency-domain  regularization  technique  is 
demonstrated  by  using  this  technique  to  penalize  those 

components  of  Hinv(£)  for  which  /<lHz  and 
/  >50  Hz  .  The  estimated  impulse  response  is  shown  in 
Figure  15.  Again  in  this  case  the  average  mean  squared 
error  was  0.1210  cm  which  is  very  close  to  the 
unregularized  value  of  0.1209  cm  which  is  again 
evidence  of  an  ill-posed  problem. 

CONCLUSION 

This  paper  presented  an  approach  to  full-vehicle 
simulator  control  which  accounts  for  nonlinearities  in 
vehicle/simulator  system.  To  improve  on  the  standard 
linear  methods  of  Cryer  et  al.  a  composite  linear  and 
nonlinear  approach  to  inverse  system  identification  was 
developed.  The  linear  method  developed  operates  in  the 
time-domain  and  employs  regularization  to  express 
preferences  for  model  behavior  in  either  the  time-  or 
frequency-domain.  The  nonlinear  portion  of  the  model  is 
designed  to  learn  the  nonlinear  deviations  from  the  linear 
model  and  the  vector-valued  support  vector  machine  is 
used  to  learn  this  deviation.  Two  approaches  are 
presented  which  describe  how  to  iteratively  improve  the 
drive  command  estimate  which  are  called  predict-correct 
and  observe-adapt.  The  linear  regularized  techniques 
are  demonstrated  for  a  2x2  system. 


Figure  13.  Plot  of  least  squares  result  of  Hinv(£) 
for  unregularized  least  squares  learning  (y= 0 ). 


Figure  14.  Plot  of  least  squares  result  of  Hmv(£) 
for  time-domain  regularization  with  settling  time 
defined  by  /u  =  0.1 ,  ks  =  400  and  y  =  l . 


Figure  15.  Plot  of  least  squares  result  of  Hinv(£) 
for  frequency-domain  regularization  with  settling 
time  defined  by  /low=lHz,  /up=50Hz  and 

7=  1- 
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DEFINITIONS,  ACRONYMS,  ABBREVIATIONS 

ASD  Auto  Spectral  Density. 

CSD  Cross  Spectral  Density. 

DFT  Discrete  Fourier  Transform. 

FFT  Fast  Fourier  Transform. 

FIR  Finite  Impulse  Response. 

FRF  Frequency  Response  Function. 

FVTR  Full  Vehicle  Test  Rig. 

HMMWV  High-Mobility  Multipurpose  Wheeled  Vehicle. 

MIMO  Multiple  Input  Multiple  Output. 

RDECOM  Research  Development  and  Engineering 
Command. 

RMS  Root  Mean  Squared. 

SISO  Single  Input  Single  Output. 

SVM  Support  Vector  Machine. 

SVR  Support  Vector  Regression. 

SYS-ID  System  Identification. 

TARDEC  Tank  Automotive  Research  Development 
and  Engineering  Command. 


